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, We present an exact calculation of the mean first-passage time to a target on the surface of a 

(N . 2D or 3D spherical domain, for a molecule alternating phases of surface diffusion on the domain 

(~| ' boundary and phases of bulk diffusion. The presented approach is based on an integral equation 

, which can be solved analytically. Numerically validated approximation schemes, which provide 

more tractable expressions of the mean first-passage time are also proposed. In the framework of 
\^ ' this minimal model of surface-mediated reactions, we show analytically that the mean reaction time 

^SJ . can be minimized as a function of the desorption rate from the surface. 



e 

I 



PACS numbers: 



INTRODUCTION 



The kinetics of many chemical reactions is influenced by the transport properties of the reactants that they involve 
[U, In fact, schematically, any chemical reaction requires first that a given reactant A meets a second reactant B. 
This first reaction step can be rephrased as a search process involving a searcher A looking for a target B. In a very 
dilute regime, exemplified by biochemical reactions in cells Q which sometimes involve only a few copies of reactants, 
the targets B are sparse and therefore hard to find in this search process language. In such reactions, the first step 
of search for reactants B is therefore a limiting factor of the global reaction kinetics. In the general aim of enhancing 
'■^ the reactivity of chemical systems, it is therefore needed to optimize the efficiency of this first step of search. 
^ I Recently, it has been shown that intermittent processes, combining slow diffusion phases with a faster transport, 
O , can significantly increase reactions rates 0, HI • A minimal model demonstrating the efficiency of this type of search, 
I ^i ' introduced to account for the fast search of target sequences on DNA by proteins [1] is as follows (see also 0413 )■ 
The pathway followed by the protein, considered as a point-like particle, is a succession of ID diffusions along the 
^-H DNA strand (called sliding phases) with diffusion coefficient Di and 3D excursions in the surrounding solution. The 
^ \ time spent by the protein on DNA during each sliding phase is assumed to follow an exponential law with dissociation 
£0 ' rate A. In this minimal model, the 3D excursions are uncorrelated in space, which means that after dissociation from 
DNA, the protein will rebind the DNA at a random position independently of its starting position. Assuming further 
. that the mean duration of such 3D excursions T2 is finite, it has been shown that the mean first-passage time at the 
target can be minimized as a function of ri = A^^, as soon as the mean time spent in bulk excursions is not too long. 
Quantitatively, this condition writes in orders of magnitude as T2 < / Di, and the minimum of the search time is 
obtained for ri ~ T2 in the large L limit. Note that in this minimal model, where the time T2 is supposed to be a 
fixed exterior parameter, bulk phases are always beneficial in the large L limit (i.e. allow one to decrease the search 
time with respect to the situation corresponding to ID diffusion only). 
. ^ In many practical situations however, the duration of the fast bulk excursions strongly depends on the geometrical 
. properties of the system (Tl| - [T^ and cannot be treated as an independent variable as assumed in the mean-field (MF) 
^ model introduced above. An important generic situation concerns the case of confined systems [lB4i3|> involving 
transport of reactive molecules both in the bulk of a confining domain and on its boundary, referred to as surface- 
mediated diffusion in what follows. This type of problems is met in situations as varied as heterogeneous catalysis 
[iMIIi'l: or reactions in porous media and in vesicular systems [H, [3[20'|. In all these examples, the duration of bulk 
excursions is controlled by the return statistics of the molecule to the confining surface, which crucially depends on 
the volume of the system. This naturally induces strong correlations between the starting and ending points of bulk 
excursions, and makes the above MF assumption of uncorrelated excursions largely inapplicable in these examples. 

At the theoretical level, the question of determining mean first-passage times in confinement has attracted a lot 
of attention in recent years for discrete random walks [2l| - |25j and continuous processes [26l - l29| . More precisely, the 
surface-mediated diffusion problem considered here generalizes the so-called narrow escape problem, which refers to 
the time needed for a simple Brownian motion in absence of surface diffusion to escape through a small window 
of an otherwise reflecting domain. This problem has been investigated both in the mathematical [2^, l30l - [3^ and 
physical [33l - l36| literature, partly due to the challenge of taking into account mixed boundary conditions. The case of 
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surface-mediated diffusion brings the additional question of minimizing the search time with respect to the time spent 
in adsorption, in the same spirit as done for intermittent processes introduced above. The answer to this question 
is a priori not clear, since the mean time spent in bulk excursions diverges for large confining domains, so that the 
condition of minimization mentioned previously cannot be taken as granted, even in the large system limit. In this 
context, first results have been obtained in [37j where, surprisingly enough, it has been found that, even for bulk 
and surface diffusion coefficients of the same order of magnitude, the reaction time can be minimized, whereas MF 
treatments (see for instance [33|) predict a monotonic behavior. 

Here, we extend the perturbative results of [37] obtained in the small target size limit. Relying on an integral 
equation approach, we provide an exact solution for the mean FPT, both for 2D an 3D spherical domains, and for 
any spherical target size. We also develop approximation schemes, numerically validated, that provide more tractable 
expressions of the mean FPT. 

II. THE MODEL 

The surface-mediated process under study is defined as follows. We consider a molecule diffusing in a spherical 
confining domain of radius R (see figure [T|), alternating phases of boundary diffusion (with diffusion coefficient Di) 
and phases of bulk diffusion (with diffusion coefficient D2). The time spent during each one-dimensional phase is 
assumed to follow an exponential law with dissociation rate A. At each desorption event, the molecule is assumed to 
be ejected at a distance a from the frontier (otherwise it is instantaneously readsorbed). Although formulated for any 
value of this parameter a smaller than R, in most physical situations of real interest a R. The target is perfectly 
absorbing and defined in 2D by the arc 9 G [— e, e], and in 3D by the region of the sphere such that 9 G [0, e] where 
9 is in this case the elevation angle in standard spherical coordinates. Note that as soon as e 7^ 0, the target can be 
reached either by surface or bulk diffusion. In what follows we calculate the mean first-passage time at the target for 
an arbitrary initial condition of the molecule. 



target 




FIG. 1: Model 



III. 2D CASE 

In this section, the confining domain is a disk of radius R and the target is defined by the arc 9 (z [— e, e]. 
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A. Basic equations 

For the process defined above, the mean first-passage time (MFPT) at the target satisfies the following backward 
equations 

^t'i{e)+\[t2{R-a,e)-h{e)] = -l for e [e,27r- e], (1) 

/ I Q 1 92 \ 
^2(^ + -^ + ;:,^)^.(.,^)^-i, (2) 

where ti stands for the the MFPT starting from the circle at a position defined on the circumference by the angle 
0, and t2 for the MFPT starting from the point (r, 9) within the disk. In these two equations, the first term of the 
Ihs accounts for the diffusion respectively on the circumference and in the bulk, while the second term of Eq. ([T]) 
describes desorption events. They have to be completed by two boundary conditions 

t2{R,e) - (3) 

^1(6*) = fore^e [0,e]U [27r-e,27r], (4) 

which describe the adsorption events and the absorbing target respectively. Eq.(l2]) is easily shown to be satisfied by 
the following Fourier series 

t2(r,6l) = ao - -7- + Va„r"cos(n6i), (5) 
41^2 — ; 

n— 1 

with unknown coefficients a„ to be determined. In particular, we aim at determining the search time (ti), defined as 
the MFPT, with an initial position uniformly distributed on the boundary of the confining domain. Taking Eq.([S]) at 
r = R, we have 

R'^ ^ on r n^ f^iW if 6' € [e, 27r - e] , 
ao h> cos(n6') = < [ ' ' , (6) 



so that 



1 r^""-' 

i?"a„ = - / ti(e) cos(ne)d9 {n > 1). 



(7) 



In what follows we will make use of the following quantities: 



= i?VV-Di, (8) 



a 



and 



xsl--, (9) 



T-^\.S^^. (10, 



A 4:D 
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As we proceed to show, two different approaches can be used to solve this problem, (i) The first approach, whose 
main results have been published in Wi\ , uses the explicit form of the Green function for the two-dimensional problem 
and relies on a small target size e expansion. We recall these perturbative results below for the sake of self-consistency 
and give details of the derivation in Appendix [X] (ii) The second approach presented next relies on an integral 
equation which can be derived for ti, and leads to an exact non-perturbative solution. 
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B. Perturbative approach 

It is shown in Appendix that the Fourier coefScients of t2 {r, 9) as defined in Eq. ([5]) satisfy an infinite hierarchy of 
hnear equations, which lead to the fohowing small e expansion: 



(11) 



Note that Eq. pDl gives in particular the first terms of the perturbative expansion of the search time (ti) defined in 
([7]) and given in ^^T]. It should be stressed that since the coefficients of e*^ of this expansion diverge with oj, in practice 
one finds that the range of applicability in e of this expansion is wider for tu small. 



C. Integral equation for ti 

In this section, we first show that the resolution of the coupled PDEs ([l][2]) amounts to solving an integral equation 
for ti only. As we proceed to show, this integral equation can be solved exactly. Writing Eq. ^ as 
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and expanding its right-hand side into a Taylor series leads to 

_ i?2 (_a)fe /afet2 



— -uj'[t2{R~a,e)-t2iR,e)], (12) 



fc=i ^ 

Substituting the Fourier representation ([5]) for t2 into this equation yields 



- ^2 y ^— ^ y a^n -l)...{n-k + cos^). (14) 



9612 £,1 \2D2 402) ^ kl 

k—l n—k 

Changing the order of summations over n and fc, using the binomial formula and the expression ([7|) for give 
92ti R"^ y f aR a'^ 

UJ 



9612 2?i V2I?2 4L> 



2 pA-K-e 

— ^(a;"- l)cos(n6')y cos(n0')^i(^')d6''- (15) 



This integro-differential equation for ti can actually easily be transformed into an integral equation for ti, by inte- 
grating successively two times, which leads to 



1 / i?2 2 / ai? a2 



+ ^ y (." - /""^os(ne')ii(^')d^', 



(16) 



or equivalently to 



^.(g) = (g - e) (2^ - e - g) + y (a;" - 1) ''"^^"^^ , '^"^^"'^ ' cos(n0')V'(^')d0', (17) 

1 ^ e 



n=l 

where 



^(^)-^, (18) 
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with T defined in Eq. ([TO]) and fl = ^. Note that Eq. (HH) holds for 9 e [e, 27r - e]. When there is no desorption 
(i.e., A = 0), only the first term in Eq. (fTS]) survives, yielding the classical result [26| 

h{9)^^i9-e){27r-e-d). (19) 

The same result is obtained for a = 0, since — 1 = (1 — a/R)'^ — 1 = 0. The limit a = is in fact equivalent to 
the limit A = because, after desorption, the particle immediately returns onto the circle (a = 0) as if it was never 
desorbed (A = 0). 

D. Exact solution 

Iterating the integral equation (|17p shows that the solution il]{9) writes for 9 e [e, 27r — e]: 

oo 

i^{9) = (6* - e)(27r - e - 6*) + ^ rf„ [cos(n6') - cos(ne)] , (20) 
with the coefficients (i„ which satisfy 

cxj oo ^ oo \ 

dn [cos{n9) - cos(ne)] = ^ ( [/„ + ^ Qn,n'dn' ) [cos(n6') - cos(ne)] , (21) 

n—l n—1 ^ n' — 1 

where we introduced 



2ir-e 

C/„ = ^^ / d0'cos(n(?')(e'-e)(27r-e-0') =4— ^ Cn, 

e 

/ N / X sin(ne) , ^ , 

^„ = (tt - e) cos(ne) + — ^ (n = l,2,...). 



and 



with 



n n' + n 

sin(27ie) ' 
2n' 



cos(ne)— ^ - cos(n t)^^ L , , ^^"(^"e) 



(22) 



1 - a;" 

Q„,„. = ^/e(n,n') (n,n' = 1,2,...), (23) 



/27r-e 
cos(n0)(cos(n'0) — cos(n'e))d6' 

/ r , cos(n'e) sin(rie) sin((n' + n)e) sin((rt' — n)e) 
= (1 - d«,n') 2 — 



(24) 



2n 

Since Eq. (PT|) should be satisfied for any 9 G [e, 27r — e], one gets d = Vl(TJ + Qd), from which 

d„ = r![(/-r!Q)-iC/]„ (n = l,2,...). (25) 

Since 

oo 

{I-^QY^ = Y,(SIQ)\ (26) 

Eq. (PO)) with the (i„ given by Eq. (I^5t can be seen as a series in powers of fl, whose n-th order coefficient is explicitly 
written in terms of the n-th power of the matrix Q. 
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Note that the first term in Eq. ([20| can also be expended in a Fourier series 

Ee„[cosM)-eosM]=^ 0, otherwise, ^''^ 

where the coefficients are obtained by multiplying this equation by cosmO and integrating from to 27r: 

e„ = {n=l,2,...). (28) 

Once the d„ determined, the search time (ii) is 



2Tr 2;t-£ 

,2^1 r ,.2r, 



(29) 



E. Approximate solution 

While the previous expression of <i is exact, it is not fully explicit, since it requires either the inversion of the 
matrix / — Q,Q or the calculation of all the powers of Q. We give here an approximation of (/ — which in turn 

provides a convenient and fully explicit representation of ti. As shown numerically (see Figs. [U [31 2] and section IVl 
for more details about numerical methods), this approximation of ti proves to be in quantitative agreement with the 
exact expression for a wide range of parameters. 

This approximation relies on the fact that, in the small target size limit e — > 0, the matrix Q is diagonal, which 
mirrors the orthogonality of the {cos(n0)}„ on [0, 2tt\. More precisely, one has from Eqs. (|23I24[) : 

Qra.n — ^ra,nQn,n ~t" ), ('^O) 

and keeping only the leading term of this expansion yields 

dn « n{l - nQr,,n)-'U„, (31) 

from which we obtain the desired approximation: 



n(7r — e) cos(ne) + sin(ne) 1 — 

' COS(^7(fc j ) 

This yields an approximation for the search time: 



m ^i0~e)i2.-e-0)+ 41) E(cos(n.) - cos(n.)) J ^—^ — + »(1 - .»)/.(n, n) ' ^''^ 



F. Variations of the search time (ti) with the desorption rate A 

In this section, we answer two important questions. When are bulk excursions favorable, meaning enabling to 
reduce the search time (with respect to the situation with no bulk excursion corresponding to A = 0)? If so, is there 
an optimal value of the desorption rate A minimizing the search time? 



1. When are bulk excursions beneficial to the search? 

This question can be investigated by studying the sign of the derivative at A = 0. The mean search time from 
Eq. (1^^ can also be written as 



R 



(*i)-^(l + A,) 



27rDi 
3i?2 



{tt - - ■ (I + XQ)-'U) 



(34) 
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FIG. 2: Comparison between three approaches for computing ti{0) in 2D: the exact solution (|20l 12511 . the approximation (I32p 
and the perturbative formula (|A15|) . with D2 — 1, a = 0.01. In the first row, the other parameters are: e — 0.1, A = 1, 
and the series are truncated to A'' = 100. On the right, the absolute error between the exact solution and the approximation 
(dashed blue curve) and between the exact solution and the perturbative formula (solid red curve). The approximation is very 
accurate indeed. In the second row, the parameters are: e = 0.1, A — 1000, and the series are truncated to N = 100 for the 
exact and approximate solutions, and to N — 1000 for the perturbative solution. One can see that the perturbative solution 
is inaccurate for large values of A, while the maximal relative error of the approximate solution is below 2%. In the third 
row, the parameters are: e = 1, A = 1, and the series are truncated to N = 100. The perturbative solution is evidently not 
applicable. In the last row, the parameters are: e = 1, A = 1000, and the series are truncated to N = 100. In this case, the 
approximate solution significantly deviates from the exact one (providing mostly negative values). The perturbative solution 
is completely invalid (not shown). 
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FIG. 3: Left: In 2D, the mean time {ti) computed through Eq. (|25l I29p with = 100 as a function of the desorption rate 
A for three values of D2: D2 = 0.5 (dot-dashed blue hne), D2 — 1 (dashed green hne), and D2 ~ 5 (solid red hne). The other 
parameters are: a — 0.1 and e = 0.01. When D2 < -D2,crit ~ 0.6348... (the first case), {ti) monotonously increases with A 
so that there is no optimal value. In two other cases, D2 > -D2,crit, and (ti) starts first to decrease with A, passes through a 
minimum (the optimal value) and monotonously increases. Symbols show the approximate mean time computed through Eq. 
(|31l I29|) . One can see that the approximation accurate enough even for large values of A. Right: The derivative ^^^^ defined 
by Eq. (|35p for the same parameters. 
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FIG. 4: In 2D, the mean time (ti) as a function of e, with D2 = 1, a — 0.01, and A = 1 (left) or A = 1000 (right). The exact 
computation through Eq. H25I I29p is compared to the approximation H31I I29p and to the perturbative approach. In all cases, 
the series are truncated to A'^ = 100. For small A (A = 1), the approximate solution is very close to the exact one, while the 
perturbative solution is relatively close for e up to 1. In turn, for large A (A = 1000), the approximate solution shows significant 
deviations for the intermediate values of e, while the perturbative solution is not applicable at all. 



where Q = ^Q^l^i V = ^"iL," ■ The derivative of (ii) with respect to A is then 



If the derivative is negative at A = 0, i.e. 



27rDi 



(?7-i+2A)/ + A2Q 



U 



(35) 



(36) 



bulk excursions are beneficial to the search. This inequality determines the critical value for the bulk diffusion 
coefficient I?2,crit (which enters through r?), above which bulk excursions are beneficial: 



24 



£'2,crit 7r(7r-e)3(2ai?-a2) 7r(7r - e)3(l - a:2) 



E 



l-a;" 



(tt — e) cos(ne) + 



sin(?ie) 



(37) 



Two comments are in order: 



9 




FIG. 5: -D2,crit as a function of e is computed from Eq. (|37|l in 2D for three values of ajR: 0.01, 0.1, and 1. When e approaches 
TT (the whole surface becomes absorbing), -D2,crit diverges (not shown). In fact, in this limit, there is no need for a bulk excursion 
because the target will be found immediately by the surface diffusion. 



(i) Interestingly, this ratio depends only on ajK and e. In the limit of e — 0, one gets 



24 



' ^ ^ n— 1 



1 - a;" 



Taking next the limit ajR — > finally yields: 



A 12C(3) 



-D2,crit TT^ 



1.4615, 



(38) 



(39) 



where C stands for the Riemann ^-function. 

(ii) The dependence of the rhs of Ea. (|57)) with e is not trivial (Fig. [5])- Indeed it can be proved to have a maximum 
with respect to e, which can be understood intuitively as follows: in the vicinity of e = 0, increasing e makes the 
constraint less stringent since the target can be reached directly from the bulk; in the opposite limit e vr, the 
constraint on D1/D2 has to tend to since the target is found immediately from the surface. Quantitatively, in the 
physical limit a — ?> 0, one finds that, as soon as D2/D1 > (£'2,crit/£'i) ~ 0.68..., bulk excursions can be beneficial. 

2. When is there an optimal value of the desorption rate X mimmizmg the search time? 



If the reaction time (ti) is a decreasing function of the desorption rate A, the bulk excursions are "too favorable", 
and the best search strategy is obtained for A 00 (purely bulk search). For the reaction time to be an optimizable 
function of A, the derivative ^^^^ has to be positive at some A. This necessary and sufficient condition remains formal 
and requires numerical analysis of Eq. (j35|) . A simple sufficient condition can be used instead by demanding that the 
search time at zero desorption rate is less than the search time at infinite desorption rate: 



(ti(A = 0)) < (ii(A^c^)). 
This writes in the physically relevant limit a <C i? (using the result of 30]): 

usm{u/2) 



(40) 



TT > V-T1-' with c e = — ^ , 
D2 37rc(e 7r%/2 Jo 



•\/cos(u) + cos(e) 



zdu. 



(41) 



Finally, combining Eqs. ([37l l41]) . the search time is found to be an optimizable function of A in the limit a ^ i? if 

2 



37rc(e) D2 7r(7r 



12 



n=l 



(tt — e) cos(ne) + 



sin(ne) 



(42) 
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Knowing that c(e) = ln(2/e) + 0{e), Eq. p2)) writes in the small e limit: 

< ^ < 1^ (43) 
31n(2/e) D2 ' ^^-^l 

which summarizes the conditions for the search time to be an optimizable function of A. This case is illustrated in 
Fig. M 

IV. 3D CASE 

In this section, the confining domain is a sphere of radius R and the target is the region on the boundary defined 
hy 6 G [0, e], where 9 is the elevation angle. 

A. Basic equations 

The 3D analogs of Eqs. (HI ^ read as 

These equations have to be completed by two boundary conditions: 

t2{R,e) = h{9), (46) 
ti{0) = for 61 e [0,e], (47) 

which respectively describe the adsorption events and express that the target is an absorbing region of the sphere. 

B. Integral equation for ti 

One can search for a solution in the following form 

2 00 

t2 ir,0)^ao- + (cos 0) , (48) 



QD2 

n—1 



where stands for the Legendre polynomial of order n. Using the orthonormality of Legendre polynomials, the 
projection of t2{R^ 0) on Pm writes 

smePrnicos6)t2iR, e)de = 2 Uo - — ) 6,n,o + -i^^hr- (49) 



Knowing that 



&D2) ' 2m + 1 



^ ' ^ 1 if 61 e [0,e], ^ ^ 



the an can be written in terms of ti{9) as 

R^ 1 
ao - ^ = 2/ smeh{9)de, 



2n + l r ^''^ 
— sine Pn{cos9)ti{9)de ifn>l. 
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Taylor expanding the rhs of 



leads to 



Using Eq. (|48]) for t2 yields 
dHi 1 dti 



502 tan 61 89 



'dW t&n6~d6 

d%_ 1 dh 
'W ta.n9~d9 



aR 



= -—^u^[t2{R- a, 6) - h{R, 9)] 



k=l 



i-a)'' (dH 



(52) 



(53) 



- "»^(" - 1) . . . (n - fc + l)i?"-^-F„(cos0). (54) 

I — 1 ■ ^^ — U 



3D2 6D2 



k—l n—k 

Changing the order of summations over n and k and using the binomial formula and Eq. ()5ip for a„ finally give 



d^ti 1 dti _ R^ 
09^ ^ tan6i 09 ^ ^ Di 



/ aR 



2 00 fTT 

— ^(a;"-l)P„(cos6')(2n + l) / sin6''P„(cos6'')ii(6'')d6'', (55) 

2 . 1 Je 



where, as in previous section, x = 1 — j^. This integro-differential equation for ti can actually easily be transformed 
into an integral equation for ti, by integrating successively two times. Indeed, multiplying first both members of Eq. 
(|55p by sin 9 and integrating between tt and 9 gives 



sin 61 t[{9) 



EL 

Di 
,2 



aR 



3D2 6D2 



(cos6l + 1) 



2 °0 ^TT 

— ^(a;"-l)(P„+i(cos0)-P„^i(cos0)) j sin0'P„(cos0')ii(^')d^', 



where we have used 



/ Pn{x)dx = - \ {1 - x')P:,{x) = -J—{P^+,{x) - P„_i(a:)). 
J n(ji +1) 2n + 1 

Dividing Eq. (|56p by sin and integrating between e and finally leads to 



(56) 



(57) 



i?2 



/ aR 



In 



sin(g/2) \ 
sin(6/2); 



2n+ 1 
n(n + 1) 



(P„(cos6') - P„(cose)) / sin 6''P„ (cos 6l')ti(e')d6'', 



where we have again used Eq. ((57)) . or equivalently to 



V'(6') = In 



1 — cos ( 



» V(x" - 1) ^" (Pn(cos6') - P„(cos e)) / sin(6'')Pn(cos6i')V'(6'')d6'' 
^ nn+1 L 



1 — cos e 
with the following definitions 

in this 3D case. 



i(n + 1) 

h{9) ^ 1 , R^-{R~af 



T . 



A 



6i:»2 



(58) 



(59) 



C. Exact solution 

Iterating the integral equation Eq. (|58l) shows that the solution iIj{9) writes for 9 e [e,7r] 

' 1 — cos ( 



V'(6l) = In 



1 — cos e 



^ dn [Pn (cos 9) - P„ (cos t)] , 



(60) 
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with the coefBcients d„ which satisfy 

oo oo ^ oo ^ 

^d„[P„(cos6')-F„(cose)] =f]^f[/„+ ^ Q„,„,d„. j[F„ (cos 61) - P„ (cose)] , (61) 

n— 1 71—1 ^ n'— 1 

where we introduced the new definitions 

(x"-l)(2n+l) J^^, . ^ /l-cos^\ (l-a;")(2n + l) ^ 

^" = ? — rn / ^'"^ ^Pn(cos6l In = 2. ^..^ ^n, 

7j[n + l) J \1 — cose/ n''{n+l)'^ 

e 

ncose\ , , P„_i(cose) 
1 + — - P„ cose +^^^^-^ n=l,2,... , 
n + I J n + 1 



(62) 



and 



with 



(l-x")(2n + l) , , , 

^ '-I,{n,n') Kn' = l,2,..., 63 

n[n + 1) 



/,(n,n')= / P„(u)(P„'(M)-P„'(cose))dM. (64) 



In Appendix [Cl we compute this integral exphcitly. 

Since Eq. ([58| should be satisfied for any 6, one gets d = il{U + Qd), from which 

d„^n[{I -nQ)'^U],^ (n = l,2,...). (65) 

As in 2D, using the series expansion of (/ — ^IQ)~^, Eq. can be seen as a series in powers of il, whose n-th 
order coefficient can be explicitly written in terms of the n-th power of the matrix Q. 
Note that the first term in Eq. (|60| can also be represented as a series 



^e4P„(cos0)-P„(cose)]=<^ ^ ^ < < ^ ^gg^ 

n^i [ 0, otherwise, 

where the coefficients e„ are obtained by multiplying this equation by P„(cos0) sin 9 and integrating from to tt: 

2n + l , 

Once the dn determined, the search time (ti) can be written as 

{h)^^ y rfe sin 0^(0) = ^^|21n(^^—— j -(1 + cose)- (68) 



D. Perturbative solution 



The first terms of a perturbative expansion with respect to e can easily be obtained from the previous exact solution. 
At leading order in e, we have 

21' 

O™,- = Q'Sn + 0(e) = - ^' K... + 0(«). 

n(n + 1) 



from which 
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One finds therefore 

C30 

ij{e) = -21ne + 21n(2sin(0/2)) - 2r2^(l 
Averaging over 0, it finally yields: 

it,) « c.^T|-21n(./2) - 1 - 217 £ 4^ 

^ n—l ^ 

This result was given in [s^] without derivation. 



„, 2n+l l-P„(cos( 



n{n + l) n(n + 1) + 2f7(l - x") 



0(6) 



(1 - 2:") 



) n(n+l) + 2fi(l-j:") 



+ 0(6) . 



(71) 



(72) 



E. Approximate solution 

As earlier for the 2D case, an approximate solution can be derived. As shown numerically (see Figs. El [7l [8] and 
section IVl for more details about numerical methods), this approximation of ti proves to be in quantitative agreement 
with the exact expression for a wide range of parameters. 

This approximation relies on the fact that, in the small target size limit e — > 0, the matrix Q is diagonal, which in 
turn mirrors the orthogonality of {P„(cos0)}„ on [0, tt]. More precisely, one has 



Qm.n — ^ra,nQn,n + ^(*^ ); 

and keeping only the leading term of this expansion yields 

dn ~ 17(1 ^Qn.Ti) Uji^ 

from which 

1 — cos 



^P{e) w In 



1 — COS e 



n{n + 1) 

The mean time (ti) is then approximated as 

ih) « Lo^riln ( 

^ ' 1 Vl-cose 



2n + 1 

X-')— 7t(P„(cos6') - P„(cose)) 



(73) 
(74) 



(1 + i^)P„(cose) + ""-l^") 
n{n + 1) + - x'^){2n + l)/e(ri, n) ' ^ ' 



1 + cos e 



n^ jl- x")(2n + 1) [(1 + ^)P„(cose) + ^-^0^] ' 
2^ n{n + l) n{n + l) + n{l- x"){2n+l)I^{n,n) 



(76) 



F. Variations of the search time (ti) with the desorption rate A 

We investigate here as in the 2D case the dependence of (ti) on A. 



1. When are bulk excursions beneficial to the search? 



The sign of ^1*^'' at A = is conveniently studied by rewriting Eq. as 



In 



1 — cos e 



1 + cos e 



X{( ■ {I + XQ)-'U) 



where Q = — Q^TjT ^^'^ ^ ~ ^TB^' "^^^ derivative of (ti) with respect to A is then 



d{ti) _ R^rj 
d\ ^ 4Z)f 



In 



1 — cos e 



1 + cos e 



(rri+2A)£+A^^ 
(/ + AQ)2 



(77) 



(78) 
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FIG. 6: Comparison between three approaches for computing ti{0) in 3D: the exact solution (|60l I65[l . the approximation (I75p 
and the perturbative formula (|7H) . with D2 — 1, a — 0.01. In the first row, the other parameters are: e = 0.1, A = 1, and the 
series are truncated to N = 100. On the right, the absolute error between the exact solution and the approximation (dashed 
blue curve) and between the exact solution and the perturbative formula (solid red curve). The approximation is very accurate 
indeed. In the second row, the parameters are: e = 0.1, A = 1000, and the series are truncated to A*' = 100. One can see that 
the perturbative solution is inaccurate for large values of A, while the maximal relative error of the approximate solution is still 
small. In the third row, the parameters are: e = 1, A = 1, and the series are truncated to A*' = 100. The perturbative solution 
is inaccurate as expected for large e. In the last row, the parameters are: e = 1, A = 1000, and the series are truncated to 
A'' = 100. In this case, the approximate solution deviates from the exact one for. The perturbative solution is negative and not 
shown. 
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FIG. 7: Left: In 3D, the mean time (ti) computed through Eq. (1651 168|) with A'' = 100 as a function of the desorption rate A 
for three values of D2- D2 ~ 1.9 (dot-dashed blue line), D2 = 4 (dashed green line), and D2 — 10 (solid red line). The other 
parameters are: a — 0.1 and e = 0.01. When D2 < -D2,crit ~ 1.9997.... (the first case), (ti) monotonously increases with A 
so that there is no optimal value. In two other cases, D2 > D2,crit, and (ti) starts first to decrease with A, passes through a 
minimum (the optimal value) and then increases. Symbols show the approximate mean time computed through Eq. (|76|) . One 
can see that the approximation accurate enough even for large values of A. Right: The derivative ^^^^ defined by Eq. (|78|) 
for the same parameters. 
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FIG. 8: In 3D, the mean time {ti) as a function of e, with D2 = I, a = 0.01, and A = 1 (left) or A = 1000 (right). The 
exact computation through Eq. (|65l I68|l is compared to the approximation (|76p and to the perturbative approach. In all cases, 
the series are truncated to TV = 100. For small A (A = 1), the approximate solution is very close to the exact one, while the 
perturbative solution is relatively close for e up to 1. In turn, for large A (A = 1000), the approximate solution shows significant 
deviations for the intermediate values of e, while the perturbative solution is not applicable at all. 



If the above derivative is negative at A = 0, i.e. 



\ \ 1 — cos e J 2 / 77 

bulk excursions are beneficial to the search. This inequality determines the critical value for the bulk diffusion 
coefficient £'2,crit (which enters through -q): 

^2,crit ^ A / 2 \ _ l + cos£ \ 2{2aR~a^) 
""{l-cosej 2 J 3i?2(^.[/) 

/ 2 \ l + cose\2(l-a;2) 

(80) 



= In 



1 — cos e 



X 



^ (l-a;")(2n + l)r, , , ^ . ^ 

/ ^ , ,,4 ' (n + 1 + ncose P„(cos£ +P„-i(cose 

n— 1 ^ ' 



16 



2.5 
2 

cT 1.5 

o 

0.5 





a = 0.01 




---a = 0.1 




- - a = 1 









0.5 



1.5 



FIG. 9: D2,crit as a function of e is computed from Eq. (|80|l in 3D for three values of a/ R: 0.01, 0.1, and 1. When e approaches 
TT (the whole surface becomes absorbing), -D2,crit diverges (not shown). In fact, in this limit, there is no need for a bulk excursion 
because the target will be found immediately by the surface diffusion. In addition, i32,crit also diverges as e — >■ because a 
point-like target cannot be detected neither by bulk excursions, nor by surface diffusion in 3D. 



In the limit of e — > 0, one gets 

D2,crit (21n(2/e)-l)(l-a;2) (1 - a;")(2n + 1) 



Di 6 \^ n2(n + l)2 



In the physically relevant limit a <^ R, one has 



^ (21n(2/.)-l) 2n + l \ ^ ^^^^ ^ 2(21n(2A)-l) ^ ^^^^^ ^g^) 



Di ' 3 \^n(n+l)^ 

There are similarities and differences between the behaviors of -D2.crit hi 2D and 3D. Figure [S] shows that £'2,crit 
from Eq. (1821) is not a monotonous function of e, with the qualitative explanation which is the same as in the 
two-dimensional case. 

In contrast to the analogous Eq. (15^1) in 2D, the rhs of Eq. diverges as e — >■ 0. This divergence reflects the fact 
that a poink-like target (e = 0), which could be found within a finite time in 2D by one-dimensional surface diffusion 
on the circle, is not detectable in 3D neither by bulk excursions, nor by surface diffusion. 

2. When is there an optimal value of the desorption rate A minimizing the search time ? 

For the reaction time to be an optimizable function of the desorption rate A, it is necessary to write an additional 
condition, requiring that the bulk excursions are not "too favorable" (otherwise, the best strategy is obtained for 
A — > cxi). A sufficient condition is given by demanding that the search time at zero desorption rate (i.e., without 
leaving the boundary) is less than the search time at infinite desorption rate 

(ti(A = 0)) < (ti(A^oo)). (83) 

(ii(A = 0)) = ^(21n(2/6) - 1) + 0(6), (84) 
which writes in the physically relevant limit a R (using the result of [38]): 

(ii(A->^)> = ^(l + eln(l/6)). (85) 



Finally, this conditions leads, for small e, to 



^>-(21n(2A)-l). (86) 
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FIG. 10: The regions of optimality for the search time in 2D (left) and 3D (right). The lower bound (solid blue line) and upper 
bound (dashed red line) are given in the limit oi a <^ R and e <C 1 by Eqs. (|43ll87p in 2D and 3D, respectively. When the ratio 
D1/D2 lies between two curves, the search time {ti) is optimizable with respect to A. Above the upper bound, surface diffusion 
is preferred (A = is the optimal solution), while below the lower bound, bulk excursions may be "too favorable" (A — j- 00 may 
give the optimal solution). We recall that the lower bound was obtained from the sufficient condition H83p meaning that the 
region below the dotted line may still be optimizable. 



Combining the two conditions 



, the search time is found to be optimizable when a <^ R and e ^ 1 if 



7(^^"(^/-)-^)<i^< 2(21n(2A)-l) - 



(87) 



V. NUMERICAL RESOLUTION 



In the previous sections, we derived the closed matrix forms ([25l [65]) for the coefficients (i„ in 2D and 3D. These 
coefficients determine the angular dependence of ti{9) through the explicit representations (|20 1 \60\i in 2D and 3D, 
respectively. Although the formulas (j251 I65|) which are based on the inversion of an infinite-dimensional matrix 
(/ — flQ) remain implicit, a numerical resolution of the problem has become straightforward. In fact, one needs 
to truncate the infinite-dimensional matrix Q and vectors U and U and to invert the truncated matrix (/ — flQ) 
numerically. 

There are six parameters that determine the function ti{9): the radius R of the disk (sphere), the diffusion coeffi- 
cients Di and D2, the desorption rate A, the size e of the absorbing region, and the distance a. From now on, we set 
the units of length and time by setting R = 1 and Di = 1. Although the distance a may take any value from to R, 
the physically interesting case corresponds to the limit of small a. As we mentioned previously, the limit a — exists 
but trivially leads to searching on the surface, without intermediate bulk excursions. In order to reveal the role of a, 
we consider several values of a: 0.001, 0.01, 0.1 and 1, the latter corresponding to the specific situation when search 
is always restarted from the center. Since the diffusion coefficient D2 enters only through the prefactor T from Eq. 
(fTO|) . its influence onto the searching time ti is easy to examine. In what follows, we take three values oi D2: 0.1, 
1 and 10. The dependence of ti on the desorption rate A and the size e is the most interesting issue which will be 
studied below. 

In the previous sections, we derived several formulas for computing ti: 

• explicit representations (POI 15(1)) with the exact expressions (^51 [55)1 for the coefficients; 

• approximations (|321 1751) which were derived by neglecting non-diagonal elements of the matrix Q; 

• perturbative formulas (jlll I7ip which are valid for small e. 

For a numerical computation of the coefficients in Eqs. ()25l I65|) . we truncate the infinite-dimensional matrix Q to 
a finite size N x N and invert the matrix (/ — fiQ). In order to check the accuracy of this scheme, we compute the 
coefficients by taking several values of from 10 to 200. For Z?2 = 1, e = 0.1, a = 0.01 and A = 1, the computed 
mean time (ti) rapidly converges to a limit. Even the computation with = 10 gives the result with four significant 
digits. Note that other sets of parameters (e.g., larger values of fi) may require larger truncation sizes. 
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VI. CONCLUSION 



To conclude, we have presented an exact calculation of the mean first-passage time to a target on the surface of a 2D 
and 3D spherical domain, for a molecule performing surface-mediated diffusion. The presented approach is based on an 
integral equation which can be solved analytically, and numerically validated approximation schemes, which provide 
more tractable expressions of the mean FPT. This minimal model of surface-mediated reactions, which explicitly 
takes into account the combination of surface and bulk diffusions, shows the importance of correlations induced by 
the coupling of the switching dynamics to the geometry of the confinement. Indeed, standard MF treatments prove 
to substantially underestimate the reaction time in this case Q, and sometimes even fail to reproduce the proper 
monotonicity In the context of interfacial systems in confinement, our results show that the reaction time can 
be minimized as a function of the desorption rate from the surface, which puts forward a general mechanism of 
enhancement and regulation of chemical reactivity. 



Appendix A: Another approach in 2D 



In this Appendix, we describe another theoretical approach which relies on the explicit form of the Green function 
of the Poisson equation in 2D case. In particular, the perturbative analysis for small e becomes easier within this 
approach. 

Considering t2 as a source term in the Poisson type equation ([1]) with absorbing conditions at 9 — e and 6 — 2ti ~ e 
whose Green function is well known 39] , ti writes 



tm = 



1 



cj sinh(2cj(7r — e)) 



27r-e 



sinh(w(0< — e)) sinh(a;(27r — e — 6'>)) 



d6i', (Al) 



and the notations 6'< = min(6', &) and 6'> = max(^?, B'). 
Injecting Eq. ([5|) into Eq. (|A1[) leads to 



A sinh(2w(7r — e)) 



1(0, 



1 + A ao 



[R-af 



A ^ a™(i?-a)"I(m,( 



where, for m integer, 
I(to,6') 

so that 



sinh(tj(0< — e)) sinh(ti;(27r ■ 



.))cos(m6i')d6l' 



uj^ + m' 



{cos{ni9) sinh(2a;(7r — e)) — 2 cos(TOe) sinh(w(7r — e)) cosh(a;(6' — tt))) 



E 



-|- 



(R-aY 2 



''^'^ = A+"°-^^ 

cosh(a;(6' - tt)) / 1_ 
cosh(w(7r — e)) I A 

Substituting Eq. (|X4l) into Eq. (O gives 



^tanh(c.(. - .)) ^ _ . + 1 - - + A f - (l - 1) 

WTT TT TT \4:D2 V TT / 4^2 



{R - a)™ cos(m6') 

OO 

y J R~arcos{me) 



m—l 



\ 2 °° 

A. ^ 



(A2) 



(A3) 



(A4) 



TT ^ — ' CJ^ + m 

m—l 



,(i?-ar^^, (A5) 
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and 



nn \ \ AD2 
XuP_ ( _ sm{{m~n)e) ^ _ sm{{ni + n )e) 

2 tanhfwfTT — e)) ^ / , , nsm(ne) \ , , 

, o ox w cos ne , , / \, , A6 

7r(a;2+n2) \^ ^ ' tanh(a;(7r - e)) / ^ ' 

S = l + \(ao- + y _ „)m cos(me). (A7) 

^ ^ m— 1 

Eq. (|A5|) can be rearranged into 

— tanh(^(7r-£)) \ _ fl f {R-af YX / e tanh(^(7r - e)) 



where 



— 7 — ^ 7T- { tanh(cj(7r — e)) cos(ne) H sin(ne) ) q;„, (A8) 

n— 1 



and Eq. into 



i?" - f 1 - -) (i? - a)" «« = 

(ao - 77^ + T) — 2 sin(ne) + — ^ tanh(w(7r - e)) cos(ne) 

— 7 — T / '^rn — 5 r- cos(TOe) (w cos(ne) tanh(tj (tt — e)) — n sin(ne) ) 

^ ' m— 1 

a; + m"^ m — n ^-^ + to^ m + n ' 



1. Particular case A = 



In this case, the previous equations can be solved exactly, leading to 



4/^2 O L*! TT 

and 

2 i?^ n(7r — e) cos(ne) + sin(ne) 1 /, -,\ 

" 7^ 3 ^11) 

TT Di n"^ R" 

We note that the particular case a = is also described by these expressions, although it does not seem to be clear 
from Eqs. (|X5 |) -(|X6 l) . 



2. Particular case a = R 



Here again, Eqs. (jA5P - (jA6p can be solved exactly, and give : 



^2 / ^ \ I _ (. _ tanh(^(7r~e)) 



TTUJ 
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and 

2/1 i?2 \ uj ^ sin(ne) + tanh(w(7r - e)) cos(ne) 1 



tanh{a;(7r— e)) J^n * 



(A13) 



TTCJ 



3. Perturbative approach 

Expanding ao and a„ in powers of e 

ao = + + + . . . and a„ = a^f) + a^^h + a(f)£2 ^ _ ^ ^^^4) 

Eqs. (jA5l) - (jA6l) lead, after lengthy calculations, to 



ao = -^+c.2rJ(2y — ^ I -7re+ (l + 2c.^y ^ ^/ ^ ) 



i?«(cj2(l - a;") +n2) 



' (1 — x™) + / \ (1 — x™) + TO^ 

1 ^ ^ / \ m— 1 ^ ^ 

{-2 + nV + ...} 



(A15) 



Appendix B: A second integral equation satified by t\ in the 2D case 

Using Eq. ([7]) for the Fourier coefficients in Eq. (jA4p leads to a second integral equation satisfied by t\ 

rr.( cosh(w(7r-6'))\ /■^''~' , , / , cosh(a;(7r-6l)) A, 

ti0=T 1 ^ (( + / tia J0,a ^ ((Je,a da, Bl 

\ coslL(w(7r — ej) / \ cosn(a;(7r — ej) / 



where 



11 ^ 
J(0, a) = — + - y (l - I) " cos(n0) cos(na). (B2) 

n— 1 



This equation is especially well adapted to local expansions of t\{ff) in the vicinity of a ~ _R, but it can also be 
rearranged into the following integral equation, useful when a R : 

hie) = y [^^_ cosh(^(^-^)) \ ^ ^ ^ ^ %^(^')sinh(a;(0<-e))sinh(a;(27r-e-e>))d0' + 



cosh(aj(7r — e)) / sinh(2cj(7r — e)) 
cosh(w(7r — 9)) 
cosh(aj(7r — e)) 



'^~\,{a) ( J{e,a) - ^^^t^^^L^Jie,a) ) da, (B3) 



where 

oo 2 

J[e,a) ^ - V ((i _ " _ l) cos(n0)cos(na). (B4) 

n=l 

Appendix C: Computation of I^{m,n) in 3D 

In this Appendix, we provide the explicit formula for the matrix I^{m,n) in 3D case. Although technical, this is 
an important result for a numerical computation because it allows one to avoid an approximate integration in Eq. 
([M)) which otherwise could be a significant source of numerical errors. The formula (IC6p for non-diagonal elements is 
somewhat elementary, while the derivation for diagonal elements seems to be original. 



Non-diagonal elements 



The Legendre polynomials satisfy 



from which 



d_ 
dx 



(i-x^)ip4x) 



+ n{n+l)Pn{x) = 0, 





J dxPn{x) 



n{n + 1) 



(n > 0). 



and 



/ , p , ,p , , [ii-x')[p^{x)p:,ix)-pux)p:„{x)]t 

J dxPm{x)P4x) = m{m+l)-n{n+l) ^ 



Since 



we find 



(1 - x^)P^{x) = -nxPn{x) + nPn-i{x) = (n + l)xPn{x) -{n+ l)P„+i(x), 



O 

j dxPn{x) 



[xPn{x) - Pn-l{x)\\ 



n+l 



(n> 0) 



and 



/■ , „ / , X [(?"-«)a;-Pm(a;)-P„(a;)+nP„_i(a;)Pm(a;)-mPm_i(a;)P„(a;)]]^ 

/ dxPm{x)Pn{x) = ^ '-^ 

J m{m + 1) — n[n + 1) 



(m 



From the above formulas, we get 

(n - m)uPm{u)Pn{u) + (to + l)P„(u)P„_i(u) - (n + 1)P„(m)P„_i(m) 



Ie{m, n) = m- 



u = cose {m^n). 



{n + l)[m(TO + 1) - n{n + 1)] 



Diagonal elements 



We denote 



O 

= JdxP^ix). 



Using the relation 



we obtain 



Pn{x) = — -xPn-l{x) - -Pn-2{x), 



Kn = — [ dxxPn-l{x)Pn{x) - I dxPn-2{x)Pn{x) . 
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The second integral is given by Eq. (jC6l) . In order to compute the first one, we consider 

b 



= 



(x) 



d_ 

dx 



{\-x')—P^{x) 



dx 



+ n{n+l)Pn{x) 



xPnix) 



d_ 

dx 



{\-X^)—Pn-l{x) 

dx 



+ {n - l)nPn-i{x) 



2n / dxxPn-i{x)Pn{x) + 



xPr,-l{x){l - x')P!^Sx) - xPn{x){l - x')P'^_^(x) 



dx{\ - x") [P,;(a:)P„_i(:E) - Pl,^^{x)P^{x)\ . 



The last integral can be written as 



J = / dx(\ - x')Pl{x){Pr,-l{x)lP,,{x))' = [(1 - X^)Pr,-l{x)P„{x)]^ 

I 

dx{Pn-l{x)/P,,{x))[~2xP^{x)+2{l - X^)P:^{x)Pn{x)] 

b 

= [{l-X^)Pn-l{x)Pnix)]l- 2 J dx[~xPn-l{x)Pnix) + (1 - x2)P; (x)P„_i (x)] . 

a 

In the last term, we substitute (1 — x^)P^(a;) to get 

6 6 

J= [{l-X^)Pn-l{x)Pn{x)]l + 2 j dxxPn^l{x)Pn{x)~2 J dx[~nxPn{x) + nPn-l{x)] Pn-l{x) . 



Bringing these results together, we get 
6 

= 2n y dxxPn^i{x)Pn{x) 



a;P„_i(x)(l - x2)P;(x) - xP„(x)(l - x^)P;^_^{x) 



+ [(l-x2)P„_i(x)P„(a;)]^ + 2 / dxxP„^i{x)Pn{x) - 2 / dx[-na;P„(a;) + nP„_i(a;)] P„_i(x) 



so that 

b 

dxxPn^l{x)Pn{x) 



-1 



4n + 2 



(Cll) 



(C12) 



(C13) 



(CM) 



x{l - x2)[P„_i(a;)P;(x) - P„(x)P;_i(x)] + (1 - x^)Pn-i{x)Pnix) 



2n + l 



We obtain 



K„ = - 



2n- 1 
2n(2n+ 1) 



x(l - a;2)[P„_i(x)P;(a;) - Pn{x)P'^_^{x)] + (1 - x2)P„_i(x)P„(a;) 

6 



2n- 1 
2n + 1 



n - 1 


'2a;P„_2(a;)P„(a;) - nP„_i(a;)P„_ 


2(x) + (n-2)P„_3(a;)P„(x) 


n 


2(2n- 


1) 



(C15) 



(C16) 
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We can further simplify this expression by using the following identities 



(n - l)P„-i(x) - (2n - 3)xPn-2{x) + (n - 2)P„_3(a;) = 0, 
nPn{x) - {2n - l)xPn-i{x) + {n- l)P„^2(a;) = 0, 
(1 - x^)P^{x) = ~nxPn{x) + nPn-i{x), 
(1 - x2)P;_i(a;) = nxPn-i{x) - nPn{x). 



We get 



2n- 1 
~2n(2n+ 1 

n- 1 
2n(2n- 1) 

2n- 1 



2n(2n + 1) 



[Pl_^{x) + P2(x) - 2xP„(x)P„_i(a;)] + (1 - x^)Pn-i{x)Pn{x) 
(2n - l)xPn{x)Pn-2{x) - nP„_i (a;)P„_2 (x) - (n - l)P„_i(a;)P„(a;) 

nx[Pl_^{x) + P^{X) - 2xPn{x)Pn-l{x)] + (1 - X^)Pn-l{^)Pn{x) 



2n-l 
2n+l 



2n- 1 



Kn-l 



1 

2n 



((2n- + l)P„(x)P„_i(x) - nx{P^_^{x) + P^{x)) 



_ [xjP^.Ax) + P^jx)) ~ 2Pn{x)Pr^-,ix)]l ^ 2n - 1 
2n+l +2n+l 

and we know that Kq = b — a. Applying this formula recursively, one finds 

P„(&)-P„(a) 



2n+l 



where 



P„(x) = x[P^{x) + 2p2_i(.T) + ... + 2p2(x) + Po(.t)] 

- 2P„(x)P„_i(x) - 2P„_i(a;)P„_2(a:) - ... - 2Pi(a;)Po(x) + x 

n 

= ^ [2{x - l)P^{x) + [Pk{x) - Pk-i{x)]'] - (x - l)P^{x) + {x- l)P^{x) + X. 

k=l 

One can check that this function satisfies the recurrent relation 

F„(x) = F„_i(x) + x[p2(x) + P^_i(x)] - 2P„(x)P„_i(x), Fo{x) = x. 

Note that P„(±l) = P„_i(±l) = ... = ±1. 
As a result, we obtain 

r/ X uPn{u) - Pn-l{u) P„(u) + 1 

I^{n,n) = —Pn(u) — ; — -\ : — , ■u = cose. 



(C17) 



+ (C18) 
2n + 1 



(C19) 



(C20) 



n + l 



2n+l ' 



(C21) 



(C22) 
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